Subdiffusion and dynamical heterogeneities in a lattice glass model 
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We study a kinetically constrained lattice glass model in which continuous local densities are 
randomly redistributed on neighbouring sites with a kinetic constraint that inhibits the process at 
high densities, and a random bias accounting for attractive or repulsive interactions. The full steady- 
state distribution can be computed exactly in any space dimension d. Dynamical heterogeneities 
are characterized by a length scale that diverges when approaching the critical density. The glassy 
dynamics of the model can be described as a reaction-diffusion process for the mobile regions. The 
motion of mobile regions is found to be subdiffusive, for a large range of parameters, due to a 
self-induced trapping mechanism. 

PACS numbers: 75.10.Nr, 02.50.-r, 64.70.Pf 



One of the most important features of glassy dynam- 
ics is its heterogeneous character, that is the coexistence 
of slowly and rapidly evolving regions, with a charac- 
teristic size significantly larger than the molecular scale. 
Such dynamical heterogeneities, which are often due to 
a jamming phenomenon (like for instance in colloids or 
in granular materials), have been observed recently both 
in experimental systems 0- HI Q and in numerical sim- 
ulations [3- 0- In order to model these effects, a 
fruitful path, which has attracted considerable attention 
in recent years, is to introduce "Kinetically Constrained 
Models" (kcms) with a very simple (usually one body) 
hamiltonian. Steric constraints are taken into account 
through kinetic rules that forbid some transitions be- 
tween microscopic states 7]. The study of these kcms 
has emphasized the role played in the relaxation process 
by rare and localized regions (often called mobility exci- 
tations or defects) that are not completely blocked by the 
constraints [a, lij LLil • These mobility excitations diffuse 
throughout the system, eventually leading to full decor- 
relation. This simple relaxation mechanism suggests a 
somewhat universal dynamical behavior, as advocated in 
flll |. Whether or not this picture applies to all glasses, 
it relies on the assumption that mobility excitations fol- 
low a purely diffusive motion, which is not justified on 
general grounds, as we shall demonstrate below. 

In this Letter, we consider a new kcm in which the 
local variable is a continuous density, rather than a dis- 
crete variable as in most studied kcms. The interest of 
the present model, compared to previously studied kcms, 
is two-fold. First, the statics of the model is non trivial, 
and can be characterized exactly. The stationary N-body 
distribution turns out to be factorizable for all values of 
a parameter that describes the (repulsive or attractive) 
interaction between particles, which has no counterpart 
in other KCMs. Second, the presence of continuous lo- 
cal densities leads to an interesting dynamical behaviour; 
the motion of mobility excitations is found to be subdif- 



fusive for a large parameter range, due to a self-induced 
trapping mechanism (i.e., not introduced by hand in the 
model). Accordingly, dynamical heterogeneities with a 
rather rich spatial structure are observed. The model 
is defined on a lattice of arbitrary dimension d with iV 
sites. In each cell centred on the lattice site i, we de- 
fine pi as the density of particles. The dynamics, aimed 
at describing density fluctuations, corresponds to a local 
redistribution of particles across the links of the lattice. 
At each time step At = tq/N (to is a microscopic time 
scale), two neighboring sites (J, k) are chosen at random, 
and pj and pk are redistributed to become p'j and p' k : 

Pj = l(Pj+Pk), p' k = 0--q){Pj+Pk), 0<q<l- (1) 

Note that the mass pj + pk is exactly conserved at each 
step. The fraction q is a random variable, chosen inde- 
pendently both in space and in time, with distribution 
tp(q) such that ip(l — q) = yj(q), that plays the role of in- 
ternal noise in the model. We want to model the fact that 
a locally dense packing is blocked unless some low density 
cell is present in its vicinity. A simple kinetic constraint 
is to allow redistribution only when pj + pk < 2pth] in the 
following, we set p t h = 1. Thus for large densities, the 
system is no longer able to reorganize locally. One can 
expect that if the average density is high, the dynamics 
slows down dramatically and exhibits glassy behaviour. 
Note that if initially < 2 for all i, the evolution rules 
forbid any density pi > 2 at later times. 

In order to find an exact solution for the sta- 
tionary state, we choose a beta distribution tp(q) = 
T(2p)/T(p) 2 [q(l - £?)]" -1 . The case fi = 1 corresponds 
to a uniform redistribution and may be thought of as 
non interacting particles (except from hard-core repul- 
sion). The case p < 1 favors q close to zero or to one, 
and can be interpreted as an effective attraction. Con- 
versely, p > 1 favors the maximal mixing value q = 1/2, 
and mimics repulsive interactions and suppressed density 
fluctuations. Note that a model similar to (but different 
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from) the present model has been studied numerically in 
the specific case p, = 1 [H]. From the Master equation 
describing the model, one sees that a non trivial form of 
detailed balance holds 0], leading to the exact station- 
ary iV-body distribution: 
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where p is the average density. The above explicit solu- 
tion is one of the central results of this letter. It shows 
that the stationary distribution is generally not uniform 
among all available states. This must be contrasted with 
the Edwards prescription, often used for generic jamming 
problems, which only holds when p = 1. Note that the 
above steady-state distribution is obtained in the long 
time limit only if, as noted above, all the initial densi- 
ties {p° } are less than 2 and at least some links initially 
satisfy + p° k < 2 (otherwise no redistribution can oc- 
cur at all). From Eq. the 'canonical' distribution 
Pcan({Pi}) j describing a subsystem with K sites, can be 
derived in the limit 1 < K < N (T^]: 



Pcan({Pi}) = TfL ]>r^(2 - Pi) e"«] 



(3) 



with (3 = —N^dln Zn / dp. The one-site distribution is 
thus given by p(p) — cp A1 ~ 1 e' 3p , with < p < 2 and 
Zcan _ c -k _ j n or( j er ^ characterize more quantita- 
tively the glassy properties of this model, we compute 
the fraction n of mobile links, defined as links (j, k) such 
that pj + pk < 2. In the 'canonical' steady state, r\ is 
computed as: r\ = J Q 2 dpi J Q 2 dp 2 p(pi) p(/0 2 ) 0(2 -pi- p 2 ) 
and can be evaluated numerically. In the limit p — > 0, 
the kinetic constraint does not play any role and r\ — > 1. 
In the more interesting limit p — > 2, the density of mobile 
links develops an essential singularity: 
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Thus the fraction of mobile links decreases very fast for 
p — * 2, but remains finite for any average density p < 2, 
suggesting that the critical density is p c = 2, as will 
be confirmed in d = 1 below. This situation is indeed 
reminiscent of what happens in the Kob- Andersen model 
[lfij in d = 2, in which the diffusion coefficient D goes to 
with the particle density p as: InZ) ~ (1 — p)^ 1 |lfi| . 

All the above analytical results concerning static quan- 
tities are valid in arbitrary dimension d. In the following, 
we present detailed numerical simulations of dynamical 
quantities in dimension d = 1, and discuss briefly the 
case d — 2, where dynamical quantities appear to be 
more complex. The relaxation properties of the model 
can be quantified by introducing on each site i a func- 
tion 4>i(t), that we choose for simplicity to be the local 
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FIG. 1: (a) Time correlation versus t/r* for~p= 1.60 to 
1.72 (/1 = 1). (fa) r* as a function of~p, for p — 0.3 (o), 1 foj, 
2 (+) and 3 fA). (c) Stretched exponential behavior of $(f) 
with p = 2, and exponent v determined from r 2 (t) (dashed), 
(d) Dynamic scaling r* versus £ — l/rj^p, p) (same symbols 
as (b)); dashed: exponent z = 1/v (z = 2 for p = 0.3). 



persistence: 4>i (*) = 1 if Pi has never changed in the time 
interval [0, t], and = otherwise. One can then 

introduce a global correlation function <i>(i) = [(<j>i(t))], 
where (. . .) and [. . .] denote averages over the sites and 
the noise respectively. Defining the characteristic decay 
time t* through $(r*) = 1/2, one can rescale the data 
by plotting Q(t) against t/r* [Fig. ^a)]. The relaxation 
time r* is plotted as a function of p for different p in 
Fig.^b). Interestingly $(t) behaves as a stretched expo- 
nential: — ln<f>(t) ~ (i/r*) 7 , as typical in glassy systems 
[Fig. |IJc)]. The exponent 7 matches perfectly the con- 
jecture 7 = v, where v is defined from the subdiffusion 
of mobility defects: r 2 (t) ~ t 2v - see below. 

Can a non trivial cooperative length scale be associ- 
ated with the increased glassiness as p — > 2~? Clearly, 
since the stationary distribution is factorized, no static 
correlation length can grow in this regime. Thus such 
a length scale can only appear in dynamical quantities, 
such as four-point correlation functions that have been 
studied recently to describe dynamical heterogeneities 
[ToL H?l Hil Il9| . In physical terms, these dynamical het- 
erogeneities can be interpreted as the appearance of 'fast' 
and 'slow' regions, with a typical size that grows as the 
glass transition is approached. More precisely, one can 
introduce on each site i a local variable 4>* = 4>i(r*). This 
allows one to define slow sites, such that <f>* = 1, whereas 
fast sites have </>* = 0. This particular choice of (j)* en- 
sures that fast and slow regions occupy equal volumes. To 
identify quantitatively the characteristic length scale £ of 
dynamical heterogeneities, we study the fourth (Binder) 
cumulant B(~p, L) of the random variable u> = L~ d ^p - <t>* , 
where the sum is over all sites of a system of size L [18j • 
This quantity measures the 'non-Gaussianity' of w; it is 
zero for large systems, L»f, and equal to a certain con- 



FIG. 2: Inset: plot of B(p, L) as a function of~p for different 
values of L (y, = 1). The curves do not cross, except pre- 
sumably for ~p = 2. Main plot: B(p~, L) plotted as a function 
of the rescaled variable 1 / (rjL) . The collapse is very good, 
showing that the characteristic length scale I is proportional 
to the inverse of the concentration r\ of mobile links. 

stant (2/3 with Binder's normalization) for L <C £■ Using 
finite size scaling arguments, one expects B(j5, L) to scale 
as a function oi£/L. We indeed find in the present model 
that all data rescale perfectly when plotted as a function 
of r/L [Fig. |2] . Thus £ is simply the average distance 1 /ry 
between mobile links. It is then natural to look for a 
scaling relation r* ~ £ z between the relaxation time and 
the cooperative length, where z is a dynamical critical 
exponent; r* is plotted as a function of £ = 1/rj on a log- 
log scale in Fig. ^d) . Interestingly, the relation between 
t* and £ is found to depend strongly on p; for p — 0.3, 
2 and 3, the data is well fitted by a power-law for £ ^ 1, 
whereas for ji = 1, some systematic curvature appears. 

As mentioned in the introduction, the high density dy- 
namics can be described in terms of reaction-diffusion of 
mobility excitations |8j . When a redistribution occurs on 
a given link, the 'state' of the link (mobile or not) cannot 
change, due to mass conservation. However, neighboring 
links can change state since the mass associated to these 
links is not conserved by this redistribution. So it is pos- 
sible to create or destroy a mobile link when it shares a 
site with another mobile link. Denoting mobile links by 
A and immobile links by 0, one can write schematically 
these processes as (A, 0) -> (A, A) and [A, A) -> (0, A), 
whereas the simple annihilation process A — * is for- 
bidden by the conservation rule. Moreover, a chain of 
two such transitions is equivalent to the motion of a 'de- 
fect' A. At high density p, i.e. at low concentration of 
mobility, branching and annihilation of mobility excita- 
tions become rare, and mobility motion is the dominant 
relaxation mechanism (for a related discussion, see [111). 

If the motion of mobility excitations was purely dif- 
fusive, as r 2 (t) ~ Dt, with a diffusion constant D that 
remains non zero as p — > 2, then r* would scale with £ 
as r* ~ £ 2 , i.e. z — 2. This scaling law is not compatible 



FIG. 3: (a) Mean-square displacement r 2 (t) for fj, = 0.3, 0.6, 
1, 1.4, 2, 3, 4 and 5 (top to bottom) and ~p — 1.75; dashed: 
slope 1. (b) Exponent v defined by r 2 (t) ~ t 2u (o), annealed 
(dotted) and quenched (dashed) barrier models predictions, 
(c) Rescaled correlation Gi(k,t) against kt v (p = 1.60, n = 
2). (d) Rescaled susceptibility X^ify/xt versus t/r* for ~p = 
1.50, 1.55, 1.58 and 1.60 (fi = 2); dots: slope 2v. 



with the data shown in Fig. Qfd), at least when fi > 1. 
Such a discrepancy could come from a critical dependence 
of the diffusion constant D on p. This mechanism is in- 
deed at play in the Fredrickson- Andersen (FA) model in 
d = 1, leading to z — 3 [Io| . However, the mechanism 
operating in the present model is different, and related to 
a genuine subdiffusive motion of individual defects. We 
show in Fig. Ota) the mean-square displacement r 2 (i) of 
mobility excitations for several values of p. The motion 
of mobility is found to be subdiffusive for a whole range of 
parameter fi, which we estimate to be \x > 1 (finite time 
effects induce some uncertainty on this threshold value). 
The exponent v characterizing the asymptotic power law 
regime r 2 (t) ~ D('p)t 21 ' is shown in Fig. Bib) where v is 
independent of p at high enough density. This relation 
suggests £ 2 ~ D(jd)t* 2u ; neglecting the dependence of D 
on p leads to z = 1/v. Using the measured values of v, 
one can compare this prediction with the direct evalua- 
tion of z on Fig. ^d) for fi = 2 and 3 (dashed lines). 
The agreement is quite good, although some discrepan- 
cies appear, presumably due to the density dependence 
of D. For ii = 0.3, one recovers the standard exponent 
2 = 2. The subdiffusive motion of mobility also accounts 
very well for the stretched exponential behavior of $(t). 
At short time (t <C r*), the correlation should behave 
as 1 - $(i) - r(t)/£, i.e. as t v . Thus $(i) is well ap- 
proximated by a stretched exponential with exponent v, 
which indeed matches perfectly the numerical data even 
for t ~ t* - see Fig. QJc). Interestingly, the dynamics 
of mobility can be mapped onto a one-dimensional bar- 
rier model [20j . since the processes (0, A, 0) — > (0, A, A) 
and (0, A, 0) — > (A, A, 0) involve waitingtimes which be- 
come broadly distributed when p — > 2 [l3j . The absence 
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FIG. 4: Direct visualization of dynamical heterogeneities with 
the variables </>*, for ~p — 1.50, 1.65, 1.75, in a system of size 
N = 100 2 (p, = 1). Only black sites (4>* = 0) have changed 
state between t = and t = r* . The typical size of both types 
of regions clearly increases with density. 

of quenched disorder suggests to consider an annealed 
model, where each random barrier is drawn anew after 
a jump. This would lead to v = l/fi for fi > 2 [3 III, 
which does not match the numerical data [Fig.[3{b)]. But 
as the environment of a mobile link is frozen for times of 
the order of t*, the quenched barrier model might be 
more appropriate. Indeed, the corresponding prediction 
v = 1/(1 + fi) for fj, > 1 is in rather good agreement with 
our numerics. Therefore, the existence of slow regions 
induces non trivial correlations that mimic the presence 
of quenched disorder, which is actually 'self-generated'. 

A more detailed description of the heterogeneities is 
provided by the following (four-point) correlation: 

G 4 (r,t) = [(Mt)&+r(t)) - (Mt))(4>i+r(t))}, (5) 

or its Fourier transform Gi(k,t). One expects the 
rescaled correlation G^k.t) = G 4 (£;, ^/(i^G^r — 0,i)) 
to scale as a function of k t v |2lJ . The corresponding data 
are plotted on Fig. GIc), showing a reasonable collapse. 
Integrating G4 (r, t) over r yields the dynamical suscepti- 
bility X4(0i which encodes some important information 
on the dynamics of the system 21]. Fig. |3fd) displays 
Xi{t)/x*4 versus t/r* for different p, where xl is the max- 
imum of X4(i), and r* has been determined from $(£). 
The resulting collapse is rather good, although very long 
times would be needed to test quality of the collapse for 
t S> t* . At short time (t <C r*), the data converges very 
slowly when p — > 2 to the predicted power law X4(i) ~ t 2l/ 
|2l|. indicating strong sub- leading corrections. 

Let us briefly discuss the qualitative behavior of the 
model in two dimensions -a fuller account will be given 
in Ref. . It is interesting to visualize the variables (j)* 
for a given realization of the dynamics. This is done on 
Fig. |U for densities p = 1.50, 1.65 and 1.75. The typical 
size of the fast and slow regions is clearly seen to increase 
with p. A strong asymmetry appears: slow regions are 
essentially compact, whereas fast regions seem to develop 
a fractal structure when their size increases. Numerical 
results (not shown) indicate that in dimension d = 2, 
the cumulant -B(p, L) cannot be simply rescaled by the 
typical distance r^ 1 ! 2 between mobility defects. Instead, 
an approximate rescaling can be obtained using £ ~ T]~ a , 
with a w 0.35. Mobility excitations exhibit a subdiffusive 



regime for /1 > 1 and short times, before crossing over to 
pure diffusion. Such a cross-over is not observed in d = 1, 
where subdiffusion seems to hold at all times. 

In summary, we have analyzed a new class of kcms 
with a conserved, continuous density field, and a param- 
eter /i accounting for interaction between particles. Al- 
though non trivial, the statics of the model can be worked 
out exactly. Our main result is that the motion of mo- 
bility excitations is subdiffusive for a large range of /i, 
due to the presence of self-induced disorder. The sub- 
diffusion exponent varies continuously with n, thus rul- 
ing out the possibility to describe the glassy dynamics in 
this model in terms of standard directed percolation. An 
analytical understanding of the subdiffusive process, for 
instance through a renormalization approach , would 
be highly desirable. Besides, experimental investigations 
of the (sub)diffusion of mobility would also be of great 
interest, and systems like sheared granular cells, where 
subdiffusive motion of tracer particles has been reported 
[23J, may be promising candidates for such studies. 
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